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ABSTRACT 

We study Compton upscattering of low-frequency photons in a converging 
flow of thermal plasma. The photons escape diffusively and electron scattering 
is the dominant source of opacity. We solve numerically and approximately 
analytically the equation of radiative transfer in the case of spherical, steady 
state accretion into black holes. Unlike previous work on this subject, we 
consider the inner boundary at a finite radius and this has a significant effect 
on the emergent spectrum. It is shown that the bulk motion of the converging 
flow is more efficient in upscattering photons than thermal Comptonization, 
provided that the electron temperature in the flow is of order a few keV or 
less. In this case, the spectrum observed at infinity consists of a soft component 
coming from those input photons which escaped after a few scatterings without 
any significant energy change and of a power law which extends to high energies 
and is made of those photons which underwent significant upscattering. The 
luminosity of the power law is relatively small compared to that of the soft 
component. The more reflective the inner boundary is, the flatter the power-law 
spectrum becomes. The spectral energy power-law index for black-hole accretion 
is always higher than 1 and it is approximately 1.5 for high accretion rates. 
This result tempts us to say that bulk motion Comptonization might be the 
mechanism behind the power-law spectra seen in black-hole X-ray sources. 

Subject headings: accretion - - black hole physics - - radiation mechanisms: 
Compton and inverse Compton — radiative transfer — stars: neutron — X-rays: 
general 
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1. INTRODUCTION 

The problem of Compton upscattering of low-frequency photons in an optically thick, 
converging flow has been studied by several researchers. Blandford and Payne were the first 
to address this problem in a series of three papers (Blandford & Payne 1981a,b, hereafter 
BP81a and BP81b respectively; Payne & Blandford 1981, hereafter PB81). 

In the first paper of the series they derived the Fokker-Planck radiative transfer 
equation which took into account photon diffusion in space and energy, while in the second 
paper they studied the acceleration of photons in a radiation-dominated, plane-parallel 
shock by using the Fokker-Planck formalism developed in the first paper. In the third 
paper they solved the Fokker-Planck radiative transfer equation in the case of steady state, 
spherically symmetric, super-critical accretion into a central black hole with the assumption 
of a power-law flow velocity v (r) oc r _/3 and neglecting thermal Comptonization. For the 
inner boundary condition they assumed adiabatic compression of photons as r — > 0. Thus, 
their flow extended from r = to infinity. They showed that all emergent spectra have a 
high-energy, power-law tail with index a = 3/(2 — (3) (for free fall /3 = 1/2 and a — 2), 
which is independent of the low-frequency source distribution. 

In a extension of the work of BP81b, Lyubarskij & Sunyaev (1982, hereafter LS82) 
studied Comptonization in a radiation-dominated shock, taking into account not only the 
bulk motion of the electrons but also their thermal motion. Also Colpi (1988, hereafter 
C88) generalized the work of PB81 by including the thermal motion of the electrons along 
with their bulk motion. However, as in PB81, the inner boundary condition was taken at 
r = and thus neglected the black-hole horizon. 

The LS82 study was further extended numerically by Riffert (1988, see also the related 
work by Bekker 1988), who found a self-consistent solution for the photon spectrum, 
the velocity and the temperature profile in the infinite shock. The spectra derived there 
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exhibit all the range of features from the thermal Comptonization dominated case to the 
bulk motion dominated one. Mastichiadis & Kylafis (1992, hereafter MK92) studied the 
upscattering of low frequency photons in an optically thick, spherical accretion onto a 
neutron star due to bulk motion only. The effects of the radiation force on the flow and 
thermal Comptonization were neglected, while the neutron-star surface was considered to 
be entirely reflective. They showed that the high-energy part of the emergent spectrum is 
a power-law with energy spectral index essentially zero. No high-energy cutoff was found 
because the Compton degradation term was not included in the radiative transfer equation. 

In the present work we give an accurate numerical and an approximate analytical 
solution of the problem of spectral formation in a converging flow, which takes into account 
the inner boundary condition, the free fall motion and the thermal motion of the electrons. 
The inner boundary is taken at a finite radius and the spherical surface there is considered 
to be anything from fully reflective to fully absorptive. The fully absorptive inner boundary 
mimics a black-hole horizon. No relativistic effects (special or general) are taken into 
account in this work and thus our spectra, though very suggestive, make instructive sense 
only, without being directly comparable with observations. However, our work derives 
for the first time the spectral-energy power-law index as a function of accretion rate. It 
demonstrates by using numerical and analytical techniques that the extended power laws 
are present in the resulting spectra. Thus, it constitutes a significant step towards an exact 
treatment of Comptonization in matter accreting onto compact objects. Some of the results 
presented here have already appeared elsewhere (Titarchuk, Mastichiadis & Kylafis 1996). 

In § 2 we give a qualitative discussion of the problem to be solved and of the solution 
that we have obtained. In § 3 we solve numerically and approximately analytically the 
radiative transfer equation and discuss the properties of the emergent spectrum. Finally, 
we summarize our work and draw conclusions in B 4. 
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2. QUALITATIVE DISCUSSION 

2.1. Accretion Flow 

Let M be the rate at which matter is accreting and let its radial inward speed be 



Vb(r) = c 



1/2 



T. 



where c is the speed of light and r s is the Scwarzschild radius. This free-fall velocity profile 
makes the tacit assumption that the radiation force on the accreting matter is insignificant, 
or equivalently, the escaping luminosity is much less than the Eddington value. The 
Thomson optical depth of the flow from some radius r to infinity is given by 



r T (r) = / drn e (r)a T = m 



1/2 

(2) 



where n e (r) is the electron number density, <jt is the Thomson cross section, m = M/Me 
and Me is the Eddington accretion rate defined by 

M E B hjL = 4 " GM '"- , (3) 

c 2 arc 

where Le is the Eddington luminosity, M is the mass of the central object, m p is the proton 

mass and G is the gravitational constant. We remark here that M can be significantly 

larger than Me, while the escaping luminosity L remains much less than Le, because the 

efficiency of converting accretion energy into luminosity is significantly less than unity for 

radially accreting black holes. Thus, our assumption in equation (1) of a free-fall velocity 

profile is justified (see also § 3.6 below). 

For super-critical accretion into black holes (rh > 1), it is evident from equation (2) 
that there should be regions in the flow from which the photons escape diffusively. Of 
qualitative importance in our discussion is the trapping radius r tr (Rees 1978; Begelman 
1979) defined by the relation 

^AA TT{rfr) s ! . (4 ) 
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No quantitative importance should be given to the trapping radius, because the mean 
number of scatterings that photons undergo depends on the photon source distribution. 

For our problem it is convenient to use not the Thomson optical depth tt but the 

effective optical depth 

t = 3 Tt\t) = 3m — , (5) 

c r 

defined such that T'(r tr ) = 1. The variable r = r'/2 will replace the radial coordinate r in 
the radiative transfer equation in § 3. 



2.2. Photon Energy Gain and Loss 

The physical picture here is the following: Consider a low energy photon that finds 
itself in the accretion flow. As the photon diffuses outward, it scatters off the inflowing 
electrons. Since this, in most cases, is an almost head-on collision between a fast moving 
electron and a low energy photon, the photon gains energy on average. The problem to be 
solved can then be described qualitatively as follows: Consider a source of photons with 
frequency t'o, or a distribution around this frequency, such that hv§ <C m e c 2 . The source of 
photons is placed anywhere in the accretion flow between the inner and outer boundaries of 
the flow and with any spatial distribution. For large optical depths in the accretion flow, 
the photons diffuse in the flow and end up either in the hole or at infinity. For small optical 
depths (say of order unity), the majority of the input photons escape after a few scatterings 
without changing their energies significantly (almost coherently scattered). However, a 
small fraction of the input photons undergo effective scatterings with significant energy 
change (Hua & Titarchuk 1995). In all cases, the diffusing photons gain energy from the 
bulk and the thermal motion of the electrons. The objective is to determine the emergent 
spectrum. 
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Since the bulk motion involves high speeds, the thermal motion of the electrons has 
little effect on the emergent spectrum if kT e <C m e c 2 (see below for a more quantitative 
statement). In such a case, the asymptotic behavior of the emergent spectrum is a 
power-law. The exact value of the power-law exponent depends on the condition at the 
inner boundary and the accretion rate. 

If the inner boundary is fully reflective (MK92), all input photons escape. Furthermore, 
the repeated reflection of the photons at the inner boundary causes them to scatter many 
times with the inflowing electrons and the emergent power law is as flat as it can be. For 
very optically thick flows, the spectral energy power-law exponent is almost zero. 

In the opposite case where the inner boundary is fully absorptive, as is the case of a 
black-hole horizon, the emergent spectrum consists only of the photons that did not reach 
the inner boundary. These photons have had fewer chances to collide with the inflowing 
electrons and thus the emergent power law is steeper than in the fully reflective case. The 
power-law exponent is different than that found by previous workers and this is entirely due 
to the finite radius of the inner boundary. 

In all cases, the power-law spectrum produced by the bulk motion of the electrons has a 
cutoff at high energies due to Compton recoil. Below we investigate the relative importance 
of the bulk (converging flow) and the thermal motion of the electrons to the mean photon 
energy change per scattering. 



There are two contributions to the mean energy gain per scattering of a photon. The 
contribution is < AE^ >, \ 
and is given by (see Appendix D) 



first contribution is < AE^ >, which is proportional to the first power of v^/c (BP81a), 



,(D - 4 



< AE [ b L> >^ E— . (6) 



/;/ 



It is worth noting that an additional factor 1/3 is present in equation (6) in the case of 



plane geometry considered by LS82 and Riffert (1988). 

The second contribution is proportional to the second power of the electron velocity, 
i.e., u 2 = v 2 h + v 2 , where Vth and Vb are the thermal and bulk velocity respectively, and it is 
given by (see Appendix D and also eq. [15] of BP81a) 



,(2) 



A(kT e + m e v 2 /3) 



< AE th , b >=< AE th > + < AE^ >S3 E v e bl ' . (7a) 

m e c z 

Thus we have that 

< AE th > 4A^ 1 

< AEi 1] > + < AEi 2) > " m eC 2{4m-i +4[r r (r)/m]V3} < 6 ' ( " j 

where 

5 = 1/rhG = 51.1 x T^m' 1 , (8) 

where = kT e /m e c 2 and Tio = fcT e /(10 keV). Hence, the bulk motion Comptonization 
dominates the thermal one if mT w < 5 1 . 

The high-energy cutoff of the spectrum can be understood in the following way: As 
it is shown in Appendix D, the fractional increase in energy of a low-energy photon in its 
collision with accreting electrons of bulk velocity v\, and temperature T e is given by 

< AE>mcr « 4m- 1 + (4/3)(r T (r)/m) 2 + ^ . (9) 

E m e c z 

At the same time, the recoil effects cause a fractional decrease in the energy of the photon 
given by (see Appendix D) 

< AE >rferr E 

— S3 . (10) 

E m e c 2 y ' 

When E <^C m e c 2 , the recoil effect is negligible, resulting in a pure power-law spectrum. At 
high energies and rh > 1 the two effects are comparable and the cutoff in the spectrum 

occurs at 

E AkT 

"■' -■4m- 1 + (4/3)(r T /m) 2 + ^|, (11) 



r~^ 



m P c z m P c z 
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or E c ps m e c 2 [4/m + (4/3) (r T /m) 2 ] when /cT e <C m e c 2 and Ty is an effective Thomson 
optical depth expected to be less than m. Because of the high efficiency of the recoil effect 
at energies of order m e c 2 such photons lose their energy completely after a few scatterings 
(much less than average). 

2.3. Emergent Spectrum 

In problems such as ours, the radiative transfer equation can be solved by the method 
of separation of variables (PB81; C88). Thus, the problem is reduced to finding the 
eigenvalue and the eigenfunction of the /c-component of the solution. In other words, the 
solution for the occupation number n(r, x) is a series of the form 



ii" 



T,x) = J2c k R k (T)N k (x) , (12) 

fe=i 

where r and x are the space and energy variables respectively, Ru{t) is the eigenfunction 

of the space operator, N k (x) is the spectrum of the /c-mode of the solution and c k are the 

expansion coefficients of the space source distribution over the set of the space eigenf unctions 

{i?fc(r)}. For the A;-component of the solution [i.e., c^Rk^N^x)], the spectrum N k (x) is 

the same at any place in the converging flow. We use this modification of the method of 

separation of variables to find an approximate analytic solution (see §3.3 and Appendix E). 

We have solved the radiative transfer equation numerically because an exact solution 
of the form (12) cannot be found for all energies and because the inner boundary condition 
depends on energy and also the energy part of the radiative transfer equation (see eq. [14] 
below) depends on the space variable (optical depth). The numerical solution has shown 
that the emergent spectrum can be well approximated by an analytic spectrum consisting 
of two parts. The first is made by those input photons which escape scattered but with 
negligible change of their energies, while the second is a power law extending to high 
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energies (see eqns. [21] and [29] below). 



3. RADIATIVE TRANSFER 



3.1. The Radiative Transfer Equation and its Solution 



Now we proceed to write down and solve the radiative transfer equation. Consider 
spherically symmetric accretion onto a compact object with rate m. The inflow bulk 
velocity of the electrons is Vb and their temperature is T e . Following BP81a (their eq. [15]), 
where the bulk motion term v\ is not neglected with respect to the thermal one as it is in 
their eq. [18]), we write for the occupation number n(r, v) 

c 



dn 
~dt 



v b • Vn = V • 



1 d 



v 2 dv 



3k 



-Vn 



1 dn_ 
-v— V • v b 



-v \n + 



3 dv 

kT e + m e v 2 /3 dn 



+ j(r,i/) , 



(13) 



m e c \ h dv 

where h is Planck's constant, n(r) = n e (r)<jT is the inverse of the scattering mean free 

path, and j(r, v) is the emissivity. Substituting for the inflow velocity v b = —Vb(r)r, where 

f is the radial unit vector, and taking into account the spherical symmetry, equation (13) 

becomes in steady state 



d 2 n 
W 1 



3\ dn 1 dn \ \ d 
T+ 2jlh = 2 X d^ ~ WbX^dx 



x 4 f h ~ l n + 



dn" 
dx 



m j 

2 KC 



(14) 



where x = hv/kT e , t = r'/2 defined by equation (5) and 

5 b - 1 = 5- l f h {r) , (15) 

where f b (r) = 1 + (v b /c) 2 /(3Q). 

The dimensionless spectral energy flux F(r,v) (PB81, C88) is written in the new 
variables as 



F(t,x) = x s 



It 
3m 



x l 2 I dn 1 dn 



x- 



dr 3 dx I ' 



(16) 
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Note that the normalization factor F = 2-k x 10 23 (/cT e /keV) 3 erg cm~ 2 s~ 1 keV _1 should be 
introduced in order to express F(r, x) in flux units. 

There are two boundary conditions which our solution must satisfy. The first is that 
the total spectral flux integrated over the outer boundary should depend only on x as 
r — > 0. Thus, 

F(t ,x) =C(x)t 2 as r o ^0, (17) 

where the subscript o means outer boundary. The second boundary condition is that we 

have a boundary with albedo A at some radius r&, or equivalently at an effective optical 

depth Tb, defined by 

_ r' h 3v b (r b ) 
Tb = 2 = 2— TTin) ■ (18) 

The reason for considering a general albedo A and not only the case A = (appropriate for 

black holes) is to be able to compare our results with those obtained in previous work. The 

net energy flux through this surface is 

F(^,,)=-x»(i^) = . (19) 

/,From the setup of the problem, t& is the largest value that the variable r can obtain. 

In a recent preprint, Psaltis and Lamb (1996) rederived the photon kinetic equation for 
Comptonization, correcting several inaccuracies in the literature. Expansion of the photon 
kinetic equation in the subrelativistic regime v b /c < 1 over the spherical harmonics leads to 
equation (14) (or eq. [15] of BP81a) for the zeroth moment of the occupation number. 



3.2. Numerical Solution 

Equation (14), satisfying the boundary conditions (17) and (19), has been solved 
numerically and the emergent spectrum has been computed. The method that we used 
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is the relaxation method for the solution of the boundary problem for elliptic partial 
differential equations (e.g., Press et al. 1992). The solution of the stationary problem is 
searched as the equilibrium solution of the approriate initial value problem. 

An interesting example of a space source distribution, mimicking the external 
illumination of a converging flow by the low-energy radiation of an accretion disk, is given 
by 

S(T) = S T 2 ex V (-r T (r)/fi), (20) 

where So is a normalization constant. An external illumination produces an exponential 
source distribution over the converging flow. The factor r 2 oc r~ 2 takes into account the 
space dilution of the radiation, while \i is the cosine of the angle of incidence at r = r s . 
The physical meaning of this space source distribution is the following: Consider a source 
of soft photons outside the converging flow and at a distance Rq *^> r s from the center 
of the black hole. The source illuminates from the outside the accretion flow and the 
radiation penetrates with an exponential probability to radius r, where it is isotropized by 
scattering. Thus, the scattered photons at radius r constitute the source of soft photons 
that get Comptonized in the accretion flow. We note here that this does not violate the 
condition of no incoming radiation at the outer boundary. This is a standard procedure for 
the reduction of the illumination problem to the problem of source distribution within the 
atmosphere (e.g., Sobolev 1975). 

EDITOR: PLACE FIGURE |TJ HERE. 

A specific example of the emerging spectra, related to the source space distribution 
(20) and to the source spectral distribution as a blackbody spectrum, is shown in Figure 1. 
Two of the emergent spectra are presented there. The first one (solid line) is a solution of 
equation (14) with the appropriate boundary conditions (17) and (19), and 5b as given by 
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equation (15). The second spectrum (dashed line) is for the case where the second order 
effect of the bulk motion has been neglected in the radiative transfer equation. Neglecting 
the (vb/c) 2 term in equation (14) transforms it into equation (18) of BP81a. It is evident 
from Figure 1 that the emergent spectral energy flux consists essentially of two components: 
The first one is due to the majority of the input photons which escape without significant 
energy change, while the second one is the flux of those relatively few photons which 
undergo effective repeated scatterings and produce an extended power law. At the same 
time, the inclusion of the second order bulk effect produces an identical shape of the hard 
tail (for energies higher than 10 keV) but with a larger photon contribution which, in 
addition, is extended to higher energies. This is because the second order effect of the bulk 
motion (vb/c) 2 increases the effectiveness of photon upscattering. 

3.3. Analytic Approximation of the Solution 

In view of the power-law that is evident in the emergent spectrum in Figure 1, it is 
natural to seek an approximate analytic solution of equation (14) where a power law is a 
part of the solution. First, a "solution" of the form (12) (so called Green's function) for 
equation (14) with the boundary conditions (17) and (19) is derived (for a justification 
see § 3.4 below). Three assumptions are made in order to separate the energy and space 
variables (i.e., to present the solution as series [12]). 

The first one is that we neglect the bulk motion term (vb/c) 2 in equation (14) which is 
a function of the space variable r. It is shown in Figure 1 that the inclusion of this term in 
the equation does not change the spectral shape of the high-energy part of the spectrum. 

The second one is that we assume that the Green's function (i.e., the spectrum for 
5— function energy injection) is a power law up to a high energy cutoff. In Appendix B the 
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energy Green's function is found (eq. [Bl]) and the validity of the power-law assumption is 
demonstrated in Appendix C (see also Fig. 4a, curve 1). 

The third one is that, in order to get the complete system of eigenfunctions, we keep 
the same boundary condition (eq. [A4] of Appendix A) for all of them. 

In Appendix A the space eigenvalue problem is solved under these assumptions and 
the transcendental equation (A6) is derived for the spectral index determination. The 
approximate analytic solution for the spectral flux can be written as the series (see eqns. 
[12] and [16]) 

00 Ci, 

F tot (0,x) = 4nr 2 F(0,x) = C F £ -^-/ fe (x,x ) , (21) 

k=l ZA k ~ <3 

where 

C F = 107rr s 2 (3/2) 3 / 2 m 5 / 2 , 

Cfc are the expansion coefficients of the space source distribution S(t) given by equation (B9), 
A| are the eigenvalues of the boundary value problem (A4-A8), while Ik(x, Xq) = x 3 Gk{x, Xq) 
is the energy Green's function (with Gk(x, Xq) the occupation number Green's function of 
the Comptonization equation [Bl]). In the case where xo <C 1 (i.e., the input photon energy 
E < kT e ), I k is given by 



h 



) , for X< Xq] 

XqJ /r,r,\ 

T(2n k )\x J Jo 
Here bk = «fc(«fc + 3 + 5), /ife = |[(<5 — 3) 2 + 8A|5] 1//2 , and the spectral index is determined by 

a k = 1 -[(5-3y + 8\l5Y/ 2 -(3 + 5)/2, (23) 

The normalization of Ik(x, Xq) is chosen in such a way as to keep the photon number equal 

to 1/xq. 

The above relations correspond to a monochromatic photon injection with energy Eq. 
In order to get the resulting spectrum for an arbitrary source spectrum g(Eo) one has to 
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convolve this with the energy Green's function. Thus, 

roo 

J k (E)= I k (E,E )g(E )dE . (24) 

Jo 

EDITOR: PLACE FIGURE g HERE. 

We have checked the numerical versus the approximate analytic solution for a wide 
range of rh and blackbody spectra of temperature T^. In Figure 2 we show a specific 
example of the emergent spectral energy flux for m — 3, kT e = 1 keV, kT^ = 1 keV and 
rb = r s . The solid curve is the emergent spectrum found from the numerical solution 
of equation (14) without the v\ term and the dashed curve is our approximate analytic 
solution. As in Figure 1, here also the emergent spectrum consists of two components: 
a) The bulk of the soft input photons that escaped from the accretion flow without any 
significant change of their energies, b) A relatively small fraction of the soft photons that is 
comptonized to high energies and thus forming the extended power law. It is worth noting 
the discrepancy between the two solutions at the high energy cutoff. Our assumption 
regarding a power law shape of the Green's function is certainly violated around the high 
energy cutoff E c , where the recoil effect becomes important. 

3.4. Properties of Emergent Spectra 

The upscattered power-law spectrum with the high energy cutoff is accurately described 
by the first term of series (21). The contribution of the different modes is clearly seen in 
Figure 3. The high energy part of the spectrum is determined by the first mode only, but 
the low-energy bump is formed as a result of the summation of all modes. 

EDITOR: PLACE FIGURE |3| HERE. 
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In a realistic situation where radial and disk accretion into a black hole occur 
simultaneously, expression (20) is only a very rough approximation, but this introduces no 
qualitative differences in the emergent spectra. A space source distribution of soft input 
photons different than that in equation (20) affects only the bump of soft photons seen and 
leaves unaffected the power law. The deeper the soft photons are injected, the smaller the 
bump. The shallower the injection is, the fewer the photons that get upscattered into the 
power law (Fig. 4a-b). 

EDITOR: PLACE FIGURE | HERE. 

In Figures 4a and 4b the analytic emergent spectra are presented for different space 
source distributions. In Figure 4a, curve 1 is for a distribution proportional to the first 
eigenfunction (eq. [A3], with (7 = 0). Curves 2 and 3 are for a distribution appropriate for 
external illumination (eq. [20]) with illumination cosines \l — 1 and \i = 0.3, respectively. It 
is evident that the spectral slope of the high-energy power law is independent of the space 
distribution of low-energy photons and, as we shall show in § 3.5, it is determined by the 
mass accretion rate rh only. In addition we want to draw the attention of the reader to curve 
1 in Fig. 4a. This curve represents the Green's function (22). It is clearly seen that a pure 
power law describes the spectrum perfectly over six energy decades. It can also be proved 
rigorously that, despite the presence of the exponential in front of the integral in equation 
(22), the power law can extend to x ^> 1. As it is shown in Appendix C, the integral in 
equation (22) is proportional to e x for x & 5, while it is proportional to x ai+3+5 for x ^ S. 
Thus, for kT e <C m e c 2 , the power law extends to energies of order m e c 2 independently of the 
electron temperature. For very high energies, x ^ S, the spectrum exhibits an exponential 
cutoff. 
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Now we proceed with the mathematical analysis of the properties of the emergent 
spectra. The source term j(r,x) in the right hand side of equation (14) can be represented 
in the factorized form j(r,x) = S(r)g(x) without any loss of generality. For example, S(r) 
can be given by equation (20) and g(x) can be a blackbody spectrum. If one neglects for the 
moment the energy change in the process of scattering (i.e., assuming coherent scattering), 
then the radiative transfer equation can be written in operator form as 

L T n coh = --[g(x)/x 3 ]S(r) , (25) 



where 



d 2 ( 3\ d 



is the space operator (see eq. [14]). The solution of this equation with the appropriate 
boundary condition (A4) is given by the series 

mg(x) ~ c k 
n coh (r) = -75—3- 2^ T2 R k{T) • (27) 

Z X fc=l A k 

Thus, the flux reads [cf eq. [21]) 

0^) = yE|^)- (28) 

1 fc=i A fc 

In the case of bulk motion dominated Comptonization, \\ ^> 3/2 for all k > 2 (see, e.g., 
Fig. 5), and as a result of this the corresponding energy Green functions Ik{x,Xo) are steep 
and act approximately as 5— functions on the input spectrum. Thus we can approximate 
the emergent upscattered spectrum as the following combination 



f W (0,x)^ Cl 



Ji(x) g(x) 



+ F t c : t h (0,x), (29) 



A? -3/2 Xj 

where J\{x) is given by the convolution shown in equation (24). Equation (29) means 
that the emergent spectrum can be written approximately as the sum of two components. 
The first component is the difference between the first terms of the full and the coherent 
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solutions (i.e., it represents the upscattered component), while the second component is the 
fully coherent solution. 

For the input distribution (20), it is straightforward to derive that the emergent 
spectral flux is 



F tot (0,x)^F Cl 



Ji(x) g(x) 



A? -3/2 A? 
where Fq is a normalization factor in flux units, 



+ c coh g(x) \ , (30) 



° c k 



Ccoh = E TT = {^I^V 12 exp(r/ 2 ) erfcfa) , (31) 

fc=l X k 

and 

v = [m/6] 1/2 //i • (32) 

We note here that the approximation of the spectrum presented in equation (29) is 
valid even if the diffusion approximation is not valid. The only requirement is that the 
contribution of the upscattering in the emergent spectrum is efficient enough to produce 
the power-law hard tail. In such a case, the first expansion coefficient C\ should be found by 
using the first eigenfunction of the kinetic equation (and not of the Fokker-Planck equation, 
as in the case of diffusion), while the coherent component F^(0,x) is the emergent 
spectrum when no energy change occurs during scattering (Titarchuk & Lyubarskij 1995). 

3.5. Power-law Spectral Index 

In order to obtain the eigenvalues of the space operator, the transcendental equation 
(A6) is solved numerically and the results are presented in Figure 5. The transcendental 
equation (A6) is factorized in the limit of fn ^> 1 (see eq. [A8]) and the roots are determined 
from the factors. The dotted line in Figure 5 is the root of the first factor of equation 
(A8). As it also can be seen, the solid line and the dotted line perfectly match each other 
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in the wide range of I < m <, 15, while the two start deviating for larger values of rh. This 
deviation shows the limits of validity and accuracy of the factorization. 



EDITOR: PLACE FIGURE | HERE. 



The asymptotic expression {rh ^> 1) for the root of the first factor of equation (A8) is 



? 3 31-A 
A, 



n 



1/2 

(33) 



2 41 + A 

In fact this asymptotic expression for Xl describes the numerical value A^ with an error of 
less than 10% (see Fig. 5). 

Having obtained the eigenvalues of the space operator, we can calculate the energy 
spectral index a,\ of the dominant term. As it can be seen from Figure 6a, this quantity 
demonstrates a weak dependence on the electron temperature for kT e < 3 keV. It is evident 
that in the wide range of rh (2 <^ rh <, 15) we have 1 < a\ < 1.5, which is different than 
what previous workers have found (PB81, C88). 

EDITOR: PLACE FIGURE | HERE. 

For kT e <C m e c 2 such that 5> 1, equation (23) can be expanded to first order in 1/5, 
giving 

a x w 2\l - 3 , (34) 

The accuracy of the above formula is better than 10% when the temperature is 1 keV and 
it becomes better than 1% in the case where the temperature is 0.01 keV. 

As long as 5 3> 1, bulk Comptonization effects dominate the thermal ones and the 
spectral index a\ of the fundamental mode is weakly dependent on the temperature of the 
electrons. 
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Combining equations (33) and (34) we find that for fn ^> 1 and 5 ^> 1 the dominant 
spectral index becomes 



31- A 

1 2 1 + A 

with an error of less than 10% for kT e = 1 keV. 



n 



1/2 

(35) 



The dependence of the spectral index oc\ on the albedo A is displayed in Figure 6b. As 
we have already pointed out, for A = the power-law index is in the range 1 < a± < 1.5 for 
fn ;> 2. For the fully reflective inner boundary case (A — 1), oc\ — ► for m ^> 1, as it was 
found by MK92. For all values of A we have taken kT e = 1 keV and r^ = r s . 

It is worth noting that the first term of equation (21) vanishes exponentially when 
Tb — > oo, since the expansion coefficient C\ goes exponentially to zero. This happens because 
the normalization Hi of the fundamental space function r 5//2 $(— A^ + 5/2, 7/2, r) (see eqns. 
[A3], [B9], [BIO]) grows exponentially when T& ^> 1 (Abramowitz & Stegun 1970). Thus, 
series (21) converges to the solution found by PB81 and C88 when T& — > oo. 

3.6. Comptonization Enhancement Factor 

A relation can be found between the luminosity of the low-frequency source and the 
luminosity emerging from the converging flow. Let us use monochromatic, low-frequency, 
dimensionless, input luminosity with unity normalization, i.e., L = / °° S(x — x )dx = 1. 
The luminosity of any energy component Ik(x,Xo) (eq. [21]) is given by 

L k = / I k (x,x )dx . (36) 

Jo 

Hereafter the index k has been dropped for simplicity. 

One can get by using equation (22) and the integration technique (e.g., Sunyaev & 
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Titarchuk 1985), 






£(a-l) \f C + a- 1~° / 

r(or(a)r(i-«)(i-4- a K- 1 



Xn x , for a > 1 ; 



(37) 
for a < 1, 



where £ = a + 4 + 5. This result is a generalization of the result of Sunyaev & Titarchuk 
(1980) for the converging inflow case and expression (37) produces a continuous transition 
through the value a — 1. 

The important conclusion which can be deduced from equation (37) is that the 
low-frequency source flux is amplified only by a factor ~ 3 (for oc\ ~ 1.5) due to 
Comptonization in the accretion flow into a black hole. Thus, the effectiveness of the bulk 
motion Comptonization is rather weak. It is easy to show that the upscattering effect is 
negligible (L/L ph 1) for the k mode with k > 2, because the spectral indices a > 2 (eq. 
[34] holds for all modes when the temperature does not exceed 1 keV). 

In §2 we assumed a free-fall velocity profile. Since the energy gain due to the bulk 
motion Comptonization is not bigger than a factor of 3, it follows that we can safely neglect 
the effects of the radiation force in our calculations if the injected photon flux in the 
converging inflow is of order of a few percent of the Eddington luminosity. 

4. DISCUSSION AND CONCLUSIONS 

We have studied Compton upscattering in the case of thermal plasma infalling radially 
into a compact object. The flow is considered to be finite (as opposed to an infinite one 
considered in the related work of PB81 and C88) by taking as inner boundary a totally 
absorptive surface to simulate spherical accretion into a black hole. The emergent spectrum 
is calculated by solving the radiative transfer equation in the diffusion approximation with 
a soft source of input photons. 
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We have considered two different expressions of the transfer equation. In the first case 
we have included the second order effect of the bulk motion (eq. [14]), while in the second 
case we have neglected it reducing the equation to the one solved in PB81 and C88 -note, 
however, the different boundary conditions. We have solved the equation in the first case 
numerically and in the second case both numerically and approximately analytically. While, 
formally speaking, the second order bulk effect should not be neglected, by comparing our 
numerical solutions for the two cases we were able to conclude that the inclusion of this 
effect does not significantly change the properties of the emergent spectrum. Thus solving 
approximately analytically the transfer equation when the second order bulk effect has 
been neglected enabled us to study in detail the properties of the emergent spectrum. This 
spectrum consists of two components, irrespective of the assumptions on the second order 
bulk effect: The first component consists of those soft input photons that escaped without 
any significant energy change, while the second is a power law that extends to high energies. 
The power law spectral index is not sensitive to whether the second order bulk effect is 
included in the transfer equation or not (see Fig.l). Thus, by using the analytic solution 
we were able to determine the spectral index as a function of mass accretion rate rh and 
we showed that this lies between 1 and 1.5 for a wide range of rh; furthermore we showed 
that the spectral index is largely independent of the electron temperature as long as this is 
below 1 keV (Fig. 6a). Therefore we have found that the finite extent of the converging flow 
has crucial effects on the emergent spectrum for moderately super-Eddington mass accretion 
rates. For extremely large mass accretion rates, our solution tends to the solution found by 
PB81 and C88. 

Our approach cannot determine accurately the exact position of the high energy cutoff. 
When the second order effect of the bulk motion is included, the numerically obtained 
spectra show that the cutoff is of order of the electron rest mass. This is further confirmed 
by Monte-Carlo calculations (Kylafis & Litchfield 1997, in preparation). If the second 
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order bulk effect is neglected, then the cutoff is of order of m e c 2 jrn as this is confirmed by 
numerical and analytical estimates. At any rate, both formalisms show qualitatively the 
same result, i.e., that bulk motion Comptonization produces a power law which extends to 
very high energies, of order the electron rest mass. 

Our spectra hint at what may be happening in candidate black-hole X-ray sources. We 
find bulk motion Comptonization as an exciting possibility for the formation of extended 
power laws in the spectra of black-hole X-ray sources. 
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A. THE SPACE BOUNDARY CONDITIONS AND THE EIGENVALUE 

AND EIGENFUNCTION PROBLEM 

In this Appendix we solve the radiative transfer equation (14) taking into account the 
boundary conditions (17) and (19), and under the assumptions stated in § 3.3. We are 
looking for a general solution of the form 

n(x,r) = R(t)N(x) . (Al) 

The space part of equation (14) can be written as 

J2 p jp 

44- (r + 3/2)^ + ^ = 0, (A2) 

dr z ar 

and its solution is (see, e.g., Abramowitz & Stegun 1970) 

R{r) = (7$(-A 2 , 3/2, r) + CY 5 / 2 $(-A 2 + 5/2, 7/2, r) , (A3) 

where $(a,b,z) is the confluent (or degenerate) hypergeometric function and C and C are 
arbitrary constants. Here A 2 are the eigenvalues of the problem, which we determine below. 
One can show that the boundary condition (17) implies that C = 0; otherwise, as r — > 0, 
F(r,x) oc r 1 / 2 , contrary to equation (17). 

In Titarchuk & Lyubarskij (1995) it is proved that a power law is the exact solution of 
the radiative transfer kinetic equation up to the exponential turnover where downscattering 
due to the recoil effect becomes important. Thus we assume that the energy part of the first 
term is a power law N(x) oc x~ £l up to the exponential cutoff. Note that the occupation 
number spectral index E\ is related to the index ot\ of the spectral energy flux by E\ = oc\ + 3 
(see, e.g., eq. [19]). The boundary condition (19) implies for R\(t) at r = T& that 

^-t*i = °> ( A4 ) 

where 

e-e 3(l-A) rr>f/° fA5) 

£i ~ £i ~ w+a) U • (A5) 
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This boundary condition, along with equation (A3), implies that the eigenvalues are 
the roots of the equation 

u ~ i Tb ) $(-A2 + 5/2 ' 7/2 ' r&) + -^ t ^~ x2 + 7 / 2 > 9 / 2 > r ») = ° • ( A6 ) 



This can be further written as 



$(-A 2 + 5/2,7/2,r 6 ) 



'5 £i \ 5-2A 2 $(-A 2 + 7/2,9/2,r fe ) 



. (A7) 



^2 3 7 7 >(-A 2 + 5/2,7/2,7^ 

Using the asymptotic form of the confluent hypergeo metric function $(o, b, z) for large 
argument z ^> 1 (Abramowitz and Stegun 1970, eq. [13.1.4]) 

$(a,b,z) = ^\e z z a - b [l + OCT 1 )], 

r(o) 

and representing T(a) as the product 



n—l 



T(a) = T(a + n)/ \J (a + i) for n = 1, 2, ..., 



j=0 



2 



we obtain the asymptotic form of equation (A7) which is valid for r^ ^> A 

(ft-x^fi (*£!_*) _„. (A8) 

^From the first term of the product we get e-y ~ 3 which implies (eq. [A5]) 
£i^3 + 1.5[(1 - A)/(l + A)]^/^) 1 / 2 . The other roots 5/2, 7/2, ..., (2n+l)/2, ... 
are related with the second, third and nth term of the above product respectively. 

Two more steps are required to calculate the emergent spectrum. The first one is the 
determination of the energy Green's function and this is presented in Appendix B. The 
second one is the proof that the spectrum is a power law up to the exponential cutoff; this 
is done in Appendix C. 
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B. THE SOLUTION OF THE ENERGY EQUATION AND THE 
EXPANSION COEFFICIENTS OF THE SPACE SOURCE 

DISTRIBUTION 

In this Appendix we obtain the solution of the energy part of the radiative transfer 
equation (14). This is derived after substituting the factorized expression (Al) for the 
occupation number into equation (14) and it reads 



1 d 

x 2 dx 



4 dN k 

x — h N k 



- 8x—^ = 7 fc A" fc - g(x) , (Bl) 

dx 



dx 
where g(x) is the spectral part of source function. The boundary conditions are 

x 3 N k — ► when x — > or x — ► oo . (B2) 

In order to find the Green's function G k (x, Xq), the spectral part of the source function g(x) 
is assumed to be proportional to a delta function, i.e., 

g(x) = S(x — x )/x 3 , (B3) 

and the differential operator of equation (Bl) should be rewritten in the self-adjoint form 
(d/dx)\p(x)(d/dx)]. The integrating factor f{x) = e x x 2 " 5 allows us to get such a form [i.e., 
p(x) = x 2 f(x)] with p(x) = e x x 4 ~ 6 . The reduction of equation (Bl) to the self-adjoint form 
is followed by integration of the left and right side of the equation from xo- to x$ + and thus 
the constant D of the symmetric form of the Green's function is obtained (see below). The 
Green's function is given by 

( DN 1 (x)N 2 (x ), Hx<x ; 
G k (x,x ) = { (B4) 

(DN 1 (x )N 2 (x), iix>x , 

where Ni(x) and N 2 (x) are the solutions of the homogeneous equation (Bl) [i.e., with 
g(x) = 0]. The solutions N\(x) and N 2 {x) satisfy the boundary conditions (B2) at x — > 
and x — > oo respectively. Hence the constant D, expressed through the Wronskian, 

W[iViOr ),A^o)] = N 1 (x )N^(x ) - N[(x )N 2 (x ) , (B5) 
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reads 

D = fM (B q\ 

x 3 p(x Q )W[N 1 (x ),N 2 (x )} ' l ; 

The solutions Ni(x) and N 2 (x) are expressed in terms of the Whittaker functions 
(Abramowitz & Stegun 1970) as 

Nl (x) = x (5 - 4 y 2 e-^ 2 M 2+s/U{s _ 3)2+4lk]1/y2 (x) , (B7) 

and 

N 2 (x) = x^ 2 e- x / 2 W 2+s/U{s _ 3r+4lk]1/V2 (x) . (B8) 

The product of p(x )W[N 1 (xo), N 2 (x )] is independent of x and is equal to 
—T(2ak + 4 + 5)/T(ak). In the limit of low energy source photons (xq <C 1), the 
Green's function G k (x, x ) is transformed into equation (22). 

The expansion coefficients Ck of the space source photon distribution S(r), where 
j/kc = S(r)g(x)/x 3 (see eq. [14]), are determined by using the orthogonality of the 
eigenfunctions (A3). Thus we obtain 



Ck 



f Tb [e -r^_xl + 5 /2 7/2, T)S(r)]dT/H k (n) , (B9) 

JO 



where 



H k {r b ) = r dre~ T r^ 2 ^ 2 {-\ 2 k + 5/2, 7/2, r) . (BIO) 

Jo 
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C. THE ASYMPTOTIC FORM OF THE ENERGY GREEN'S FUNCTION 

In order to study the asymptotic form of the solution (21) when kT e <C m e c 2 , the 
integral of equation (22) should be calculated for indices u; = 5 + 3 + ai^>lby the steepest 
descend method. Introducing a new integration variable t' — to, the integral 



J(x, a,u) 



is rewritten in the form 



e-t'ix + tTt'^dt' 



dt 



f 00 dt 

J(x, a, uo) = a a exp [ — <jc(i)l — • 
Jo ' t 



(CI) 



(C2) 



Here and below the subscript 1 is omitted for simplicity. The exponent of the integrand in 
equation (C2), 



u 



has a maximum at 



(p(t) = a[t \n(x + at) — hit] 

a 



u + a — x +[{uj + a — x) 2 + 4xa] 



1/2 



The second derivative of ip(t) is 



<p"(t) = a 



2a 



uoa 



(x + at) 2 t 2 _ 
Finally the integral J(x, a, uj) is given by the steepest descend formula 



(C3) 



(C4) 



(C5) 



J(x, a,u) = a' 



t exp[-ip(t ) 
to 



2k 



V"(*c 



1/2 



(C6) 



This formula provides very high accuracy for any index a > 1, but for a < 1 the integral 
should be transformed by using integration by parts as follows 



J(x, a, uj) = [J(x, a + 1, uj) — u>J(x, a + l,u — l)]/a 



(C7) 



; From the above we find that 



J(x, a, u) = e x T(uo + a) = e x T(2fj,) , for x < u , 



(C8) 
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and 

J(x, a, u) = r(a)a; w+3+a , for x > u . (C9) 

Thus, the spectral energy power law E~ a extends up to a high energy cutoff. Beyond 
this energy, the spectrum falls off exponentially (or has a hump) of the the form 

C h x s+3 exp(-E/kT e ). 



-30- 



D. ENERGY GAIN AND LOSS PER SCATTERING 

There is a standard procedure to estimate the energy gain and loss per scatering 
(e.g., Prasad et al. 1988). In order to make these estimates one has to find the response 
of the energy operator of the kinetic equation (see eq. [14]) to a delta-function injection 
of photons. The dimensionless photon energy density in this case is x 3 n = x5(x — xn), 
where x = E/kT e , E is the photon energy, T e is the temperature of the electrons and n 
is the occupation number. The energy change per scattering Ax can then be found by 
multiplification of the kinetic equation by x 3 and integration over x. Thus we have 



Ax 



1 * dn , nt d 

■—x — + x® fa- 
in ox ox 



w^+a 



dx 



where /& = 1 + (t>&/c) 2 /(30), and 6 = kT e /m e c 2 . Integration by parts gives 



(Dl) 



Ax 



-+4[0 + K/c) 2 /3] 



m 



x 3 ndx — / x A ndx . 
o Jo 



Substituting x 3 n = x5(x — x ) in this equation we obtain 



Ax 



- + 4[0 + (v b /c) 2 /3] 



E n 



m 



m P c z 



x i 



(D2) 



(D3) 



where E is the photon energy before scattering and Xn = E /kT e . The first term on the 
right hand side of equation (D3) is the Comptonization to first order in the electron bulk 
velocity, the second term is to second order in the electron velocity (thermal and bulk) and 
the third term gives the photon energy loss due to the recoil of the electron. 
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E. SEPARATION OF VARIABLES 

We remind the reader how the method of separation of variables works in the case of a 
general linear operator equation of the form 

L x n + L T n = S(r)g(x), for < t < 17, and < x < oo, (El) 

with boundary conditions l n = at r = and l Tf n = at r = t&, where r and x are the 
space and energy variables respectively, and the boundary operators l , l Tb are independent 
of the energy variable (in the case considered in our present paper, l n depends on energy). 
The eigenfunction R k {r) of the space operator (L T R k + \\R = 0, l Rk = l Tb Rk = 0) always 
exists for the very wide class of operators considered in equation (El). The eigenf unctions 
form a complete set, i.e., the space part of the source function S(t) can be expanded in a 
series S(t) = TjC k R k (r), with c^ the expansion coefficients. Having in mind this expansion 
and also that n(r, x) is looked for as a series 

n(r,x) = J2c k Rk(r)N k (x) , (E2) 

fc=i 

where Nj,(x) is the spectrum of the k-mode of the solution, one can see that the whole 

problem is reduced to the solution of the equation for the k-energy component N k (x). 

Namely, 

Zc k R k (r)[L x N k - \\N k - g{x)\ = 0, (E3) 

from which one obtains 

L x N k -\ 2 k N k -g(x) = 0, (E4) 

because of completeness of the {R k (T)}. 

In general, when the internal boundary operator l x T depends on energy (as it is in our 
case, see eqs. [16, 19]), the above method of separation of variables cannot be applied. But 
in the particular case of BP81a (their eq. [18]; see also our eq. [14] with 5b = S), a power 
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law is a solution in the range of dimensionless energies x$ < x < 5 (x C 1, 1 <C <J). This 
can be checked by simple substitution of the product R(r)x~ a into equation (14) (where 
5b has to be replaced by 5) and into the boundary conditions (eqns. [17] and [19]). After 
substitution we get a system of equations for the determination of the spectral index a. 
In addition, we assume that the photon source spectrum g(x) has characteristic energies 
much less than m e c 2 (for example, for a S- function energy injection, g(E) = S(E — E ) and 
Eq <C m e c 2 ), and therefore there are no source photons with E ^> Eq. Thus, instead of 
equation (E4), we get in this case 

L x x~ a - X 2 x- a = (A 2 - \ 2 )x- a = 0, 

or 

\l-\ 2 = 0, (E5) 

for x > x and where A^ is the eigenvalue of the operator L x . The internal boundary 
operator l x T is a linear differential operator with respect to the variables r and x. The power 
law x~ a is an eigenfunction of this operator, namely l* b x~ a = l^ h x~ a . Thus 

l x Tb x- a R{r) = [llR{r)]x- a = at r = r b , 

or 

l° b R(r) =0 at r = r b , (E6) 

and 

l R(r) =0 at r = 0. (E7) 

As it is proved in Appendix C, a power law is an eigenfunction of the operator L x for 
energies up to a high energy cutoff. Combining equation L T R + X 2 R = and equations 
(E5)-(E7) we get a complete system of equations for the determination of the spectral index 
a, the eigenvalue A 2 and the eigenfunction. Practically, the correctness of this approach 
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can be clearly seen from Figure 2, where the numerical and the analytical solutions are 
compared. 

It is not our goal to go into the deep mathematical details about the uniqueness of this 
solution. We can only remind the reader that there is a general theorem for the existence 
and uniqueness of the solution for the boundary problem of the partial differential equation 
of the second order. This theorem is valid for a wide class of differental operators L x , L T , 
l x I 
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Fig. 1. — Plot of the emergent spectral energy flux versus photon energy for rh — 3, kT e = 1 
keV, kT^ = 1 keV, A — 0, \x — 0.3 and the radiative transfer equation is solved numerically. 
The dashed line is for the case where the second order bulk effect is excluded (see BP81a, 
their eq. [15]), while the solid line includes this effect (eq. [14] of the present paper). 

Fig. 2. — Plot of the emergent spectral energy flux versus photon energy. The solid line is 
identical to the dashed one in Figure 1, while the dashed line corresponds to the approximate 
analytic solution. 

Fig. 3. — Plot of the emergent spectral energy flux versus photon energy. Here kT e = 2.5 
keV, A = 0, Tg/Tbi, — 3, rh — 3. The emergent spectrum is calculated as a convolution of 
series (21) with a blackbody source distribution. The space source distribution is given by 
equation (20) with So = 1 and \x = 0.3. The solid line is the emergent spectrum. The dashed 
lines represent the contributions to the spectrum by the first nine terms of series. 

Fig. 4. — a) Plot of the emergent spectral energy flux versus photon energy for a 
monochromatic injection at energy E = 10~ A kT e and for three different space source 
distributions. Curve 1 is for a space source distribution proportional to the first 
eigenfunction. Curves 2 and 3 are for the space source distribution given by equation (20) 
with \x — 1 and \x = 0.3 respectively. The other parameters are the same as in Figure 3. 
b) Same as in a), but for a blackbody input spectrum with kT^ = 0.833 keV. 

Fig. 5. — Plot of the first four eigenvalues as functions of the accretion rate rh in the case 
of A = , kT e = 1 keV and r& = r s . The solid curves correspond to Af, A|, A| and A4 from 
bottom to top. The dashed curve gives the approximate analytic value of A^. 

Fig. 6. — a) Plot of the energy spectral index ol\ versus rh for r^ = r s and kT e = 0.01 keV 

(solid line), kT e = 1 keV (dotted line) and kT e = 3 keV (dashed line). 

b) Plot of the energy spectral index a± versus rh for different values of the albedo A. Here 
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kT e = 1 keV, rb = r s and A equal to (solid line), 0.5 (dotted line) and 1 (dashed line). 
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